Antimatter spectra from a time-dependent modeling of supernova remnants 
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We calculate the energy spectra of cosmic rays (CR) and their secondaries produced in a super- 
nova remnant (SNR), taking into account the time-dependence of the SNR shock. We model the 
trajectories of charged particles as a random walk with a prescribed diffusion coefficient, accelerating 
the particles at each shock crossing. Secondary production by CRs colliding with gas is included as 
a Monte Carlo process. We find that SNRs produce less antimatter than suggested previously: The 
positron/electron ratio J- e + / T e + +e - and the antiproton/proton ratio Tp/J-p+p are a few percent 
and few x!0~ J , respectively. Both ratios do not rise with energy. 



Introduction — Measurements of the antimatter frac- 
tion of cosmic rays (CR) provide not only insight into 
CR physics itself [l|, as e.g. their propagation in the 
galaxy, but are also valuable probes for cosmology and 
particle physics. In particular, the annihilation of dark 
matter (DM) leads to an equal injection rate of matter 
and antimatter particles into the Galaxy, while the CR 
flux from astrophysical sources is matter-dominated. A 
possible way to detect DM is therefore to estimate care- 
fully the expected antimatter fluxes from astrophysical 
sources and to search then for any excess @. 

The PAMELA collaboration presented recently results 
of their measurement of the positron fraction in CRs, 
which is rising rapidly from 10 to 100 GeV |3|- At the 
same time, the antiproton ratio measured by PAMELA 
declines above 10 GeV Q, consistent with expectations. 
The conventional estimate for antimatter fluxes from as- 
trophysical sources uses as only production mechanism of 
antimatter the scattering of CRs on interstellar gas [l|. 
As discussed e.g. in Ref. [3], the energy dependence of the 
Galactic diffusion coefficient, D oc E s with 5 = 0.5 — 0.6, 
is inconsistent with an increase of the antimatter fraction 
with energy. By contrast, the spectral shape of fragmen- 
tation functions leads quite naturally to such a rise in the 
case of DM annihilations. 

The DM interpretation of the PAMELA excess faces 
however several difficulties Q: First, the required rate of 
positron production is larger than expected for a stable 
thermal relic. As a consequence, either the annihilation 
rate has to be enhanced, or the DM particle should be 
unstable with the appropriate life-time. Second, in gauge 
boson or quark fragmentation, positron, antiproton and 
photon production are tied together and thus one has 
to postulate a DM particle annihilating only into elec- 
trons and muons. More importantly, assuming antimat- 
ter production by diffusing CRs as the only astrophysical 
source for antimatter falls short: Since electrons lose fast 
energy, the high-energy part of the e~ + e + spectrum 
should be dominated by local sources as nearby pulsars, 
as pointed out already 20 years ago [6]. Moreover, elec- 
tromagnetic pair cascades in pulsars result naturally in 
a large positron fraction together with a "standard" an- 



tiproton flux. 

More recently, supernova remnants (SNR) were put 
forward as an alternative astrophysical explanation for 
a rising positron fraction Q: Positrons created as sec- 
ondaries of hadronic interactions in the shock vicinity 
participate in the acceleration process and, according to 
Ref. p| , should thus have a flatter energy spectrum than 
primary electrons. It was estimated that the resulting 
positron fraction can explain the PAMELA excess and 
rise up to 50% at higher energies @, while subsequently 
a similar mechanism for antiprotons was suggested in 
Ref. [8|. Since shock acceleration in SNR is expected to 
be the main source for Galactic CRs [9] , such a scenario 
has also important consequences for the interpretations 
of CR data as, e.g., the boron-to-carbon ratio [10j |. 

The present work examines the production of sec- 
ondary p and e + in SNRs, improving on previous stud- 
ies [HI m two respects: First, we use a Monte Carlo 
(MC) approach calculating the trajectory of each parti- 
cle individually in a random walk picture. This makes 
it easy to include interactions and the production of sec- 
ondaries. Second, our approach allows us to include the 
time (and spatial) dependence of relevant parameters de- 
scribing the evolution of a SNR as, e.g., the shock radius 
and its velocity, the magnetic field or the CR injection 
rate and to test their influence on the CR spectra. We 
should also stress what are not the aims of the present 
work: We do address neither the problem of acceleration 
from a microscopic point of view nor consider any feed- 
back of CRs on the shock or the magnetic field. Although 
the latter processes are important to obtain accurate CR 
escape fluxes, we shall show that our simplified treatment 
leads to an upper limit on the secondary fluxes. 

Simulation procedure — Shocks around SNRs are sup- 
posed to be collisionless, with charged particle scatter- 
ing mainly on inhomogeneities of the turbulent magnetic 
field. We model such trajectories by a random walk in 
three dimensions with step size lo(E) determined by an 
energy- dependent diffusion coefficient D. Diffusion close 
to the shock is usually assumed to proceed in the Bohm 
regime with the mean free path Iq proportional to the 
Larmor radius R L . Thus D(E) = d /3 = c 2 p/(3ef B B), 
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where Jb denotes the ratio of the energy density in the 
turbulent and in the total magnetic field. We neglect the 
coupling between CRs and the turbulent magnetic field, 
assuming that a layer with Bohm diffusion extends far 
enough into the up-stream region. For a constant mag- 
netic field, CRs do not escape but are confined in the 
SNR, corresponding to an "age-limited" scenario for the 
CR flux from SNRs. 

We describe the evolution of the shock in the rest frame 
of the SNR. Then the (yet unshocked) up-stream region 
is at rest, V\ = 0, and has the density of the surround- 
ing interstellar medium (ISM), p\ = prSM- Assuming a 
strong shock with Mach number M. ^> 1, the shocked 
down-stream region flows with the velocity V2 — 3v s h/R 
and has the density p2 = Rpi ■ Here, R denotes the com- 
pression ratio R — (7 + 1)/(7 — 1) = 4 for a mono-atomic 
gas with 7 = 5/3. We account for the flow, adding in 
the down-stream region on top of the random walk an 
ordered movement of the particle with velocity V2 that 
is directed radially outwards. Thus a particle trajectory 
evolves during the time step At — Iq/c as 

x(t + At) = x(t) + v 2 At <d(r ah - r) e r + 1 , (1) 

where lo denotes a random step, r s h the position of the 
shock, and i9(x) the step function. 

Crossing the shock front, particles are accelerated. We 
neglect that the relative energy gain £ = (E^+i — Ek)/Ek 
per cycle fc depends on the angle of the trajectory to the 
shock front, and use for simplicity that on average for a 
non-relativistic shock £ = ^(t>2 — vi) = Vsh/c. For the 
position r s h and the velocity u s h of the SNR shock we 
use the n = case of the analytical solutions derived in 
Ref. [l2[. These solutions connect smoothly the ejecta- 
dominated phase with free expansion r s h oc t and the 
Sedov- Taylor stage r s h oc i 2 / 5 . The acceleration of CRs 
is assumed to cease after the transition to the radiative 
phase at the time t max . 

As the injected particles diffuse, electrons lose en- 
ergy via synchrotron radiation and inverse Compton 
scattering, while protons can scatter on gas of the 
ISM producing secondaries that include antiprotons and 
positrons. Cross sections and the final state of pp- 
interactions are simulated using QGSJET-II [l3j], while 
we use SIBYLL 2.1 14] for decays of unstable particles. 

The last ingredient for our simulation procedure is an 
injection model. To ease the comparison with the results 
of [7], we fix the electron/proton ratio K ep at injection to 
K e p = 7x 10 -3 . As injection energy we use Eo = 10 GeV. 
In the first model used, the injection rate 

N oc r* h v? h 5(E-E )5(r-r sh ) (2) 

is proportional to the volume swept out per time by the 
shock, i.e. a = 1 (thermal leakage model [HI). I n the 
second model, the injection rate N is proportional to the 
CR pressure [16[ and a — 3. In the case of model 2, the 
fraction of particles injected very early is significantly 
larger than in model 1. 
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FIG. 1: Proton spectra (black) as function of energy for two 
different injection models and JbB = 1/iG. Additionally the 
contribution of protons staying until £ max in the up-stream 
(red) and of protons in down-stream region (blue) are shown. 



We use the following parameters to describe the SNR: 
We choose the injected mass as M e j = 4M Q , the me- 
chanical explosion energy as E snl = 5 x 10 51 erg, and the 
density of the ISM as tcism = 2 cm" 3 . The end of the 
Sedov- Taylor phase follows then as i max = 13.000 yr [l2j]. 
For the magnetic field we use B = 1/iG and /s = 1 
to ease the comparison with the stationary approach of 
Refs. [H. 

Numerical results — In Fig. [U we show the energy de- 
pendence of the proton spectra in model 1 and 2. Addi- 
tionally to the total spectra, the contribution of protons 
staying at i m ax in the up-stream region is shown in red, 
while the spectra of protons advected down-stream are 
shown in blue. The spectra of electrons are not shown, 
since they have the same shape as the proton spectra 
apart from a somewhat lower cutoff energy. While the to- 
tal energy spectra at low energies agree well with a l/£ 2 
power-law, changing the injection model leads to large 
differences at high energies. The strong dependence of 
the spectra close to E max on the injection model is ex- 
pected, since this part is sensitive to how many particles 
are injected early, when shock acceleration is most effec- 
tive. 

We switch now to the produced secondaries and show 
in Fig. [2] for injection model 2 and JbB = lpG the an- 
tiproton flux split into a part produced in the accelera- 
tion zone (A) and a part produced in the inner part of the 
SNR (B). More exactly, we define the contribution A as 
all secondaries that crossed at least once the shock. This 
contribution increases fast, since the time i acc primary 
protons stay in the acceleration zone and can interact in- 
creases as i acc oc D(E) oc E. Hence for the relative rise 
of A not the acceleration of secondaries but of primary 
protons is important. For example, the component A for 
neutral secondaries as e.g. photons, defined formally as 



3 



10" 



TTTI 1 — I I I I I 



T — I I I I II 




A+B 



1 1 ii 



10 10 10 11 



10 13 10 14 



FIG. 2: The total flux of antiprotons together with the con- 
tribution A and B in model 2 as function of energy. 



10" 



TTTI 1 — I I I I I 

model 1 
model 2 



~l — Mill 




E(eV) 



FIG. 3: The positron ratio J- e +/J- e + +e - (blue) and antipro- 
ton ratio Tp/J-p+ p (red) in model 1 (dotted) and 2 (solid). 



all the particles produced up-stream, rises with energy in 
the same way as the one of antiprotons. Since the inelas- 
ticity, i.e. the energy fraction (zp) w 0.02 transferred to 
all antiprotons is practically constant in the relevant en- 
ergy range, E /(zp) > 10 11 eV, its does not influence the 
shape of the antiproton flux [l8j]. At E h rj 2 x 10 12 eV, 
the increase of contribution A stops, the total flux retains 
its approximate E~ 2 slope and stays small in contrast to 
the result of Ref . [8] . 

How can we understand this behavior and the maximal 
value of Tp/Tp+p? We may assume in a gedankenexperi- 
ment that in each pp —> p + X interaction the most ener- 
getic antiproton carries away all the energy, Ep w (z)E p 
with (z) « 1. Then interactions just convert part of the 
p into a p flux. But since p and p diffuse and are accel- 
erated in the same way, the total p flux is not affected if 
the p or the parent proton is accelerated. Hence the to- 



tal flux of antiprotons produced in the acceleration zone 
and inside the SNR, i.e. the sum of A and B, should be 
simply the proton flux scaled down by a constant fac- 
tor. In particular, the secondary flux of the species i 
is bounded by the proton interaction depth r and the 
(spectrally averaged) energy fraction (zi) transferred to 
i. The maximal conversion rate during the life-time of a 
SNR is with cri ne i = 30 mb as inelastic pp cross section 
at 100 GeV given by r = c £ max R tt-ism find ~ 3 x 10~ 3 . 
The mean energy fraction of antiprotons (plus antincu- 
trons) is {zp) sa 0.02, so we may expect a maximal ratio 
of Fp/Fp+p ~ (zp) t ~6x 10 -5 . The obtained Fp/Tp+p 
ratio shown in Fig. [2] is indeed close to this estimate. 

The relative size of the partial contributions A and B 
can be understood considering the relation between the 
time i acc spent by protons in A, their final energy E oc 
t acc and thus the interaction depth ta in A as function 
of energy, ta oc t acc oc E. In particular, it takes all the 
life-time t max of the SNR to accelerate protons to the 
highest energies, cf. the up-stream component in Fig. 1. 
For the component B, the optical depth tb of the parent 
proton is tb oc (i max — t &cc ), which explains why the two 
components A and B sum up to a flat spectrum. Note 
that the relative normalization of component A and B 
in the stationary approach of Refs. @,H, H3 has to be 
imposed by hand, since B is formally infinite. 

The same discussion applies to the case of positrons, 
with the sole exception that the primary electron flux is 
scaled down by the factor K ep and that the energy frac- 
tion transferred to positrons is (z e +) « 0.05. The results 
of our simulation are shown in Fig. O confirming with 
the maximal value of J- e +/J- e + +e - < few % this sim- 
ple picture. Note that while above ~ 100 GeV the ratio 
J 7 e + /J 7 e + +e - from SNR starts to be larger than the con- 
ventional prediction using only secondary production on 
the ISM, it cannot explain the rise to J- e + / J- e + +e - « 10% 
at 10 GeV in the PAMELA data 0]. 

Up to now, we have discussed only our numerical re- 
sults for constant fsB = 1/iG and one may wonder if 
a "better" choice of parameters can increase the anti- 
matter fluxes. In particular, the analytical formula of 
Ref. 0, [1, EH seem to imply that the contribution A in- 
creases for weaker diffusion, i.e. larger D. However, the 
term D jv\ regulating the importance of A limits also via 
t acc oc D/v\ the maximal proton ener gy. Using a con- 
stant value Jb — 1/20 as in Ref. @, M. Ill| thus reduces 
£ max by the same factor. In the stationary approach, 
however, E max is an external parameter and by increas- 
ing £" max relative to its natural value given by i acc = i max 
one enlarges the relative contribution of A. This approach 
has been justified as a method to account in an effective 
way for amplification and damping of the magnetic field. 

In our time-dependent approach, we test this sugges- 
tion considering a time-dependent magnetic field: assum- 
ing that non-linear effects amplify magnetic fields in the 
early phase [13, with fsB = 100//G before the transi- 
tion to the Sedov- Taylor phase at t* — 240 yr, while in the 
late stage magnetic fields are damped, JbB = 1/20/iG 
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FIG. 4: Spectra of cosmic rays protons and positrons (scaled 
up by a factor 1000) together with the partial contributions 
A and B in model 2 for a time-dependent diffusion coefficient. 



at t > t*. In Fig. HJ we show for this case the pro- 
ton and positron spectra using the injection model 2. 
Protons that were injected early are accelerated up to 
few x 10 15 eV, while the bulk of CRs injected when the 
tubulcnt magnetic field is damped has a cutoff around 
10 12 eV. The contribution A to the positron flux satu- 
rates at E ~ few x 10 11 eV, i.e. at the energy expected for 
JbB = 1/20/iG. In contrast to Fig. [21 the component B 
dominates now the high-energy end of the positron flux: 
Most CRs escape from the acceleration zone during the 
transition fsB/nG — 100 — > 1/20 and those advected 
downward contribute the new high-energy extension of 
component B, while those escaping up-stream contribute 
to A but are not longer accelerated. 

Finally, we stress that the splitting between contribu- 



tion A and B is artificial and depends as well as the total 
flux on the definition of the escape flux: If the diffusion 
coefficient drops above a certain energy and/or outside a 
sufficiently small radius r s h+<5r to a value close to the one 
typical for the Galaxy, then CRs can escape up-stream 
instead of being confined down-stream. Clearly, this ef- 
fect reduces the contribution B. On the other hand, the 
bounds T e + / T e + +e - < few % and Tp/Tp +P < 6 x 10~ D 
will become stronger, since also the time for interactions 
in the acceleration zone will be shortened. Since our max- 
imal values of J-p/J-p+ p and J- e + / J- e + +e - depend only on 
t max , which is lower in escape- limited than in age-limited 
models, we conclude that the contribution of SNR to the 
observed antimatter in CRs does not lead to rising anti- 
matter fractions and is smaller than estimated in earlier 
works. 

Summary — We calculated the energy spectra of CRs 
and their secondaries produced in a supernova remnant 
using a simple random walk picture. In contrast to a pre- 
vious prediction that the positron fraction J- e + / J- e + +e - 
can rise up to 40%~50% for K ep = 7 x 10~ 3 , we found 
that the ratio levels off at a few percent. This value cor- 
responds to the expectation combining the interaction 
depth t R3 3 x 10~ 3 of a proton during the life-time of 
a SNR with the energy fraction z ~ 0.05 transferred to 
positrons. Similarly, the antiproton ratio TpjTp+p does 
not rise beyond fewx 10~ 5 . Our results suggest that anti- 
matter production in SNRs cannot explain the rise of the 
positron fraction observed by PAMELA. Since a rising 
antiproton fraction is neither expected from CR interac- 
tions with the ISM nor from pulsars, such a measurement 
could be used as a signature for dark matter. 
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